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ABSTRACT 



The CLEAN deconvolution algorithm has well-known limitations due to the restriction of locating point source model components 
on a discretized grid. In this letter, we demonstrate that these limitations are even more pronounced when applying CLEAN in the 
case of Rotation Measure (RM) synthesis imaging. We suggest a modification that uses Maximum Likelihood estimation to adjust the 
CLEAN-derived sky model. We demonstrate through the use of mock one-dimensional RM synthesis observations that this technique 
shows significant improvement over standard CLEAN and gives results that are independent of the chosen image pixelization. We 
suggest using this simple modification to CLEAN in upcoming polarization sensitive sky surveys. 
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CLEAN (H6g boml[l97l is a deconvolution algorithm that is 
widely used throughout astronomy. In CLEAN, a deconvolved 
image is constructed iteratively by locating the maximum in the 
image, adding a delta function at that location with some fraction 
of the peak brightness to a model, and then subtracting a point 
spread function (PSF) that has been shifted to the location of 
the model component and scaled by its strength. This is repeated 
until a user-defined stop condition has been reached. 

The algorithm makes an implicit assumption that the im- 
age is well described by a collection of statistically indepen- 
dent point sources. In some cases this is an adequate assumption 
and experience has shown that even in circumstances where this 
assumption is not ideal, CLEAN is still able to produce quali- 
tatively reasonable results. Because of its flexibility, simplicity, 
and speed, CLEAN has become the standard tool for image re- 
construction in radio astronomical aperture synthesis imaging. 

A limitation of CLEAN arises due to the discretization of 
the image space, and therefore one can only create model point 
sources on a regularly spaced grid. To properly remove the PSF 
pattern associated with a strong point source, it is crucial for 
the delta function representing the source in the CLEAN model 
to have an accurate position. While the PSF limits the mini- 
mum spatial scale of structures that can be resolved, the accu- 
racy with which one can measure the location of a point source 
depends on the signal-to-noise ratio. In the high signal to noise 
regime the uncertainty in the location of a source can be much 
smaller than the PSF, and to measure the source location ac- 
curately with a single model point source one would have to 
create an image with a large number of pixels within the PSF. 
Such over-resolving is at best computationally wasteful and in 
some cases may make the imaging problem completely unfeasi- 
ble. For example, even with a few pixels pe r resolution element , 
three-dimensional (3D) Faraday synthesis dBell & EnBlin 2012I) 
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images will easily exceed the available memory on a single com- 
puter. 

CLEAN is able to partially overcome this issue by describing 
a single point source by a cluster of CLEAN components around 
the true source location. The strongest model point source will 
lie on the pixel closest to the source location. Weaker model 
points will be added to neighboring grid points such that when 
the collection of model sources is convolved by a Gaussian 
restoring beam that approximates the idealized PSF, the resulting 
image will be closer to the true source position and flux than that 
of a single model point located on the grid. Nevertheless, inaccu- 
racies in image reconstructions due to the pixelization of the sky 
are a wel l known limiting factor to dynamic range when using 
CLEAN dBriggs & CornwelJl992t|p^vll999l:ICotto"n & Usonl 
l2TjollYatawattall2010l) . 

Rotation Measure (RM) synthesis dBrentiens & de Bruvnl 

2005) is a promising new technique for studying magnetic fields 
with the new generation of broadband radio telescopes. RM 
synthesis allows for the separation of multiple sources of po- 
larized emission along a line of sight when each is Faraday 
rotated by different amounts. It produces an estimate of the 
Faraday spectrum, or the polarized emission as a function 
of the Faraday depth, which is a quantity that measures the 
amount of Faraday rotation. The Faraday depth, <p, is propor- 
tional to J n e B z dz where n e is the density of thermal electrons 
and B z is the component of the magnetic field along the line 
of sight. RM synthesis is similar to aperture synthesis imag- 
ing due to the Fourier relationship between the Faraday spec- 
trum and the polarized intensity as a function of the squared 
wavelength. CLEAN has been proposed as a de convolution 
techn ique to be applied to RM synthesis imaging ( Heal deTaTI 
2009), and is often referred to as RMCLEAN in this con- 
text. Several al ternative image reconstruction techniques have 
been proposed dFricket al.ll2010t iLi et alJl20lTt lAndrecut etail 
I2012I) . but RMCLEAN has thus far been used most often 
as the deconvolution method of choice dFeain et al.l [2009: 
lHarv ev-Smith et al teolOt iMaoet al.ll2010t iBernardi et all2010t 
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In this letter we show that the limitations of CLEAN due to 
pixelization of the sky are particularly pronounced in the case 
of RM synthesis, and propose to use a Maximum Likelihood 
(ML) estimation procedure to improve the model obtained us- 
ing the standard RMCLEAN algorithm. This method is similar 
to those proposed previo usly in the context of high-fidelity aper- 



ture synthesis i maging (El-B eherv & MacPhidll980t lYata watta 
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Bern ardi et al.ll2.Qllh . It is a lso similar to the method pro- 



posed bylO'Sullivan^^I] d2012l) . but is more easily scaled to 
higher dimensions and includes fewer assumptions about the 
Faraday spectrum. Our aim is to maintain or improve accuracy 
of the standard RMCLEAN algorithm while keeping the number 
of pixels to a minimum. 

2. The algorithm 

Here we consider the case of ID RM synthesis where one is at- 
tempting to reconstruct the complex valued Faraday spectrum, s, 
from some polarized sky brightness data, d, that has been taken 
at many frequencies. We note that the following description is 
trivially extended t o 3D as would be required in the case of 
Faraday synthesis dBell & EnBlin 20121) . which combines aper- 
ture and RM synthesis into a single procedure that provides sig- 
nificantly better results than when they are performed separately. 

The Faraday depth coordinate axis (the image coordinate) 
will be denoted as <p, and the data coordinate is A 2 . For RM syn- 
thesis, the following measurement equation applies: 



d(A z ) = S(A z ) I d<K#2 2 ' ,<M +n(A l ) 



0) 



Here n is a Gaussian additive noise term and S represents a sam- 
pling function defining the discrete A 2 values at which measure- 
ments are obtained. 

Using the CLEAN algorithm, one generates a list of model 
points M - nti, </>,,• with flux values and positions given by m, and 
(pi, respectively, where i is an index over the list of Nm model 
points. The representation of the model in data space is 



(2) 



where j is an index over the Nj values of A 2 for which measure- 
ments have been made. 

We assume that the probability of measuring the full data set, 
d, given the CLEAN model (i.e. the likelihood), is a product of 
the probabilities of measuring each individual data point, dj. In 
this way, we assume the measurements to be independent of one 
another, thus making the likelihood 



P(d\M) = Y\ P{dj\M) = \\ 




(3) 



Here, as mentioned previously, we assume Gaussian distributed 
noise, and cr 2 is the variance of the noise in channel /'. 

Our approach is to adjust the CLEAN model by maximizing 
the likelihood with respect to m, and 0,-. For simplicity, we opt to 
work with the negative log-likelihood 



H(d\M) = -hP(d\M)= \Y, 



— \2 



[dj - dj) 



(4) 
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Fig. 1. A comparison of the (a) flux and (b) position errors as a 
function of image pixel size in reconstructions of a single point 
source using RMCLEAN (black line) and the ML algorithm de- 
scribed in this paper (red line). The reported errors are the aver- 
age absolute errors from 100 simulated data sets. 



The maximum likelihood solution is obtained by minimizing 
this function. There are many possible approaches to minimiz- 
ing Eq. [4] The iterative procedure that we use is presented in 
Appendix lAl 

Minimization of Eq. |4] provides an estimate for the optimal 
model parameters under the assumption that the number of point 
sources in the initial CLEAN model is appropriate. This, how- 
ever, is almost certainly not the case. CLEAN will naturally 
over-estimate the number of point sources that are required to 
model the data because it needs a cluster of sources in order to 
model a single source at an arbitrary location. 

To find the most appropriate number of free parameters in 
the model that are supported by the data, we modify H to add 
a term that penalizes additional degrees of freedom. We u se the 
so-called Bayesian Information Criterion dSchwarzll 19781 BIC) 
to suit this role. The BIC, C, is given by 



C = 2H(d\m) + 3N m lnN d . 
The algorithm proceeds as follows: 



(5) 



Start with a list of model point sources generated by 
RMCLEAN. Condense the model such that there is only one 
entry per pixel location. 

Throw away any model entries with a flux below a user- 
defined threshold. This step removes model points that result 
from cleaning too deeply or that only function to slightly 
shift a single point source location. We use a threshold of 
twice the CLEAN image noise level. This step is not neces- 
sary, but can speed up computation in case there are many 
excess point sources in the initial model. 
Take one or more steps to iteratively minimize Eq. |4] for 
the current model using e.g. the prescription outlined in 
AppendixlAl 

Attempt to merge nearby model points. Pairwise combine 
CLEAN components into a single component having a lo- 
cation at the flux-weighted mean location and a flux that is 
either the sum of the two model point fluxes or that is solved 
for using e.g. Eq. IA.7I Accept the merged CLEAN compo- 
nent if C is reduced, otherwise revert to the previous model. 
Iterate steps [3] and |4] until the fractional change in H from 
one step to the next is below a user-defined threshold. 
Obtain a residual image by subtracting the new model from 
the data. 

Convolve the new model with an idealized PSF and add it to 
the residual image. 
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Fig. 2. As Fig.Q~|but with the average errors of two sources that 
are separated by 2xFWHM (top row) and 6xFWHM (bottom 
row). 



3. Demonstration 

To demonstrate the effectiveness of the approach outlined above, 
we present a series of ID RM synthesis mock observations of 
simple single and double point-source sky models. We generate 
simulated data according to Eq. Q] including Gaussian random 
white noise and sample using a frequency coverage with 800 
channels equally spaced between 1 and 4.2 GHz. This frequency 
range is similar to that of the combination of the upgraded Very 
Large Array L- and S-band receivers and corresponds to a PSF 
with a full width at half maximum (F WHM) of th e main peak 
of ~ 40 rad/m 2 . A Clark style CLEAN (IClarkll 1980h is then per- 
formed and the resulting model is used as input for the ML al- 
gorithm. The output images from both procedures are stored for 
comparison. 

In Fig. Q] we compare the performance of standard 
RMCLEAN to that of the ML method at reconstructing the posi- 
tion and flux of a single point source as a function of pixel size. 
The source is located at roughly the same <p location at each reso- 
lution, but shifted slightly to ensure that it is always located one- 
third of the way between two pixels. The source flux is roughly 
100 times the noise level. The plotted flux and position errors are 
the average of the magnitude of those from 100 trials with dif- 
ferent noise realizations. Source position and flux are measured 
by locating the maximum of the modulus of the Faraday spec- 
trum and fitting a Gaussian to the image around this point. We 
find that for RMCLEAN the position and flux errors increase 
significantly with pixel size while those for the ML approach 
remain nearly constant. Six or more pixels per FWHM of the 
PSF are required before the positional accuracy for RMCLEAN 
matches that of the ML approach. To obtain accurate fluxes with 
RMCLEAN, 40 or more pixels per FWHM are required. 

For a sky model with a single source, RMCLEAN is found to 
work well as long as the image is sufficiently over-resolved. With 
more than one source, however, over-resolution is no longer 
effective, as shown in Fig. [2] Here we show positional and 
flux errors in the case where there are two sources along the 
line of sight. In one case, they are separated by 80 rad/m 2 
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Fig. 3. As Fig. [2] but with fixed resolution and varying source 
separation, which is given as a multiple of the FWHM of the 
PSF (40 rad/m 2 ). 



(roughly 2xFWHM of the PSF) and in the other by 240 rad/m 2 
(6xFWHM). The modulus and phase of the complex valued flux 
is the same for both sources. The errors that are plotted are the 
average of those from both sources (the errors for the individ- 
ual sources were qualitatively similar). We find that increasing 
the resolution does not reduce errors in the RMCLEAN recon- 
struction and that on average these errors are two or more or- 
ders of magnitude larger than in the maximum likelihood re- 
constructions. This behavior is due to the interactions between 
the highly structured, complex valued PSFs associated with each 
point source. The exact dependence on the error as a function of 
source separation depends on the specific sampling function, and 
therefore PSF, being used. 

We also find that the errors are larger for both approaches 
when sources are close together relative to the FWHM of the 
PSF. To test this further, we kept the image resolution fixed (at 
5 rad/m 2 per pixel) and varied the separation between the two 
sources (again, each having the same flux and phase). In Fig. [3] 
we see that below a source separation of 2xFWHM the errors for 
both methods increase dramatically. Nevertheless, at separations 
of 1 xFWHM or larger, the maximum likelihood approach fares 
significantly better than RMCLEAN. 

Even though the errors increase dramatically as source sep- 
aration decreases for both methods, we find that the ML recon- 
structions are still a significant improvement over RMCLEAN. 
In Fig. |4]we show reconstructed Faraday spectra from 500 data 
realizations of the same sky model. Although the ML results 
show a systematic reduction in flux and a shift in position rel- 
ative to the sky model, two distinct sources are clearly identi- 
fiable. This is not the case in the RMCLEAN image. We per- 
formed similar tests for hundreds of combinations of source sep- 
aration (between 1 and 6xFWHM), relative fluxes (flux ratios 
between 1 and 10), relative phases, and noise levels. In all cases 
the ML approach was an improvement over RMCLEAN, and it 
only showed significant errors in a few cases like the one shown 
in Fig. E] 



4. Conclusions 

We have shown that the well-known limitations of the CLEAN 
algorithm that arise due to discretization of the sky are particu- 
larly pronounced in the case of RM synthesis. We find that accu- 
racy of the results obtained by CLEAN depend strongly on the 
choice of pixelization, and significant over-resolution is required 
to obtain accurate measurements of flux and Faraday depth in 
the simplest case where there is a single source along the line of 
sight. Furthermore, RMCLEAN is unable to accurately recon- 
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Fig. 4. Modulus of the Faraday spectra of two point sources having similar fluxes and the same phase and that are separated by 
lxFWHM (40 rad/m 2 ). The dashed line shows the sky model after convolution with the idealized PSF. The (a) RMCLEAN and (b) 
ML reconstructions from 500 trials are plotted in red/yellow, with the color scale indicating the log of the number of reconstructions 
that pass through each location on the figure. 



struct fluxes and locations of the individual sources even with 
over-resolution when there are multiple sources of emission. 

We propose a simple algorithm to adjust the CLEAN model 
parameters using Maximum Likelihood estimation together with 
a prescription for reducing the degrees of freedom in the model 
to the minimum number that are supported by the data. We show 
that this algorithm improves upon the results obtained using 
RMCLEAN alone and provides highly accurate reconstructions 
independent of the choice of pixelization of Faraday depth space. 

Both algorithms struggle in cases where two point sources 
are located within lxFWHM of the PSF, but this merely re- 
flects the fundamental limitation of the observations to resolve 
structures on these scales. However, although errors increase as 
source separation decreases, the ML algorithm still provides sig- 
nificant improvement over RMCLEAN. 

CLEAN is ideally suited for the case where the image to 
be reconstructed is well-described by a set of independent point 
sources, but as we have shown, does not perform well even in 
these ideal circumstances when applied in the case of RM syn- 
thesis. The algorithm described herein provides significant im- 
provements over RMCLEAN alone and is easy to append to 
existing RM synthesis imaging pipelines. Therefore, we recom- 
mend its use in upcoming polarization surveys that plan to in- 
clude RM synthesis imaging. 

The ML procedure still assumes that the image is described 
by a set of independent point sources, and will likely not be ideal 
in case large-scale, diffuse emission is present. A point-source 
based image reconstruction method is nevertheless an important 
tool to have on hand. It can be used as a method for separating 
point sources from a diffuse background, which may be better re- 
constructed using other methods that do not handle point sources 
well. It may also be the optimal deconvolution method for use 
with instruments that are not sensitive to scales much large r than 
the PSF, as is the case with e.g. LOFAR dBeck et al.ll2012l) . 



Acknowledgements. This research was performed in the framework of the DFG 
Forschergruppe 1254 Magnetisation of Interstellar and Intergalactic Media: The 
Prospects of Low-Frequency Radio Observations. We thank Henrik Junklewitz 
and Marco Selig for many helpful discussions. 



References 

Andrecut, M., Stil, J. M., & Taylor, A. R. 2012, Astronomical Journal , 143, 33 
Beck, R., Frick, P., Stepanov, R., & Sokoloff, D. 2012, Astronomy & 

Astrophysics , 543, A113 
Bell, M. R. & Enfllin, T. A. 2012, Astronomy & Astrophysics , 540, A80 
Bernardi, G., de Bruyn, A. G., Harker, G., et al. 2010, Astronomy & 

Astrophysics , 522, A67 
Bemardi, G, Mitchell, D. A., Ord, S. M., et al. 2011, Monthly Notices of the 

Royal Astronomical Society ,413,411 
Brentjens, M. A. 2011, Astronomy & Astrophysics , 526, A9 
Brentjens, M. A. & de Bruyn, A. G. 2005, Astronomy and Astrophysics, 441, 

1217 

Briggs, D. S. & Cornwell, T. J. 1992, in Astronomical Society of the Pacific 
Conference Series, Vol. 25, Astronomical Data Analysis Software and 
Systems I, ed. D. M. Worrall, C. Biemesderfer, & J. Barnes, 170 
Clark, B. G. 1980, Astronomy & Astrophysics , 89, 377 
Cotton, W. D. & Uson, J. M. 2008, Astronomy & Astrophysics , 490, 455 
El-Behery, I. & MacPhie, R. 1980, Antennas and Propagation, IEEE 

Transactions on, 28, 234 
Feain, I. J., Ekers, R. D., Murphy, T, et al. 2009, Astrophysical Journal , 707, 
114 

Frick, P., Sokoloff, D., Stepanov, R., & Beck, R. 2010, Monthly Notices of the 

Royal Astronomical Society, 401, L24 
Harvey-Smith, L., Gaensler, B. M., Kothes, R., et al. 2010, Astrophysical Journal 

,712,1157 

Heald, G., Braun, R., & Edmonds, R. 2009, Astronomy and Astrophysics, 503, 
409 

Heald, G. H. 2012, Astrophysical Journal Letters , 754, L35 

Hogbom, J. A. 1974, Astronomy & Astrophysics Supplement, 15, 417 

Iacobelli, M., Haverkorn, M., & Katgert, P. 2012, ArXiv e-prints 

Li, E, Brown, S., Cornwell, T. J., & de Hoog, F. 2011, Astronomy & 

Astrophysics , 531, A 126 
Mao, S. A., Gaensler, B. M., Haverkorn, M., et al. 2010, Astrophysical Journal , 

714,1170 

Mao, S. A., McClure-Griffiths, N. M., Gaensler, B. M., et al. 2012a, 

Astrophysical Journal , 755, 21 
Mao, S. A., McClure-Griffiths, N. M., Gaensler, B. M., et al. 2012b, 

Astrophysical Journal , 759, 25 
O'Sullivan, S. P., Brown, S., Robishaw, T., et al. 2012, Monthly Notices of the 

Royal Astronomical Society ,421, 3300 
Perley, R. A. 1999, in Astronomical Society of the Pacific Conference Series, 

Vol. 180, Synthesis Imaging in Radio Astronomy II, ed. G. B. Taylor, C. L. 

Carilli, & R. A. Perley, 275 
Schwarz, G. 1978, The annals of statistics, 6, 461 

Van Eck, C. L., Brown, J. C, Stil, J. M., et al. 201 1, Astrophysical Journal , 728, 
97 

van Weeren, R. J., Rottgering, H. J. A., Interna, H. T, et al. 2012, Astronomy & 

Astrophysics , 546, A 124 
Yatawatta, S. 2010, arXiv: 1008. 1892 



4 



M.R. Bell et al.: Improved RMCLEAN reconstructions, Online Material p 1 



Appendix A: Iterative maximization approach 

We assume that the optimal model point position is shifted by some small distance, §<pj, relative to the grid point location found 
using CLEAN, <pj 

(pi =Ji + 5<pi. (A.l) 
The maximum likelihood solution is obtained by minimizing H in Eq. |4]with respect to 5<pi, which gives 



dH{d\M) 
ddcpi 



iAjtrii < m i - I dj — d'j I e 



2iA 2 4>, 







(A.2) 



where d'j is the same as dj, Eq. [2] except that the sum is over k + i. The operator % selects the real part of its argument. This 
expression cannot be solved analytically for 5<pi, but we can use it to search for the solutions iteratively. To do so we Taylor expand 
the exponential term in Eq. IA.2l to first order, thereby making H second order in §<pi. This gives 



86th o* 4 



A^niii lm* — I dj - d'j I e 2,/l ^' - 2iA *■ I d'j - dj 



= 0, 



(A.3) 



which we can rewrite in the form 



where 



A = BS<p h 



(A.4) 



i 



dj-d'j] e 2a ^' 



A^niji 



—2Ajnii I dj — d'j] e 



(A.5) 



We must also solve for an updated flux. We can again extremize Eq.|4] this time with respect to m,. We find 



(A.6) 



and thus 



i N d ~ 



7=1 



(A.7) 



In our example implementation, we solve for 5<p and m iteratively until convergence is achieved. We also attempt to merge nearby 
model components to reduce the degrees of freedom in the model according to the prescription described in Sec. [2] 

We tried other iterative schemes, e.g. solving for position and flux changes by inverting the Hessian matrix of H, but found that 
the approach given here is the most stable. 



